function [ W, H ] = NMF( duomenys, R, niter, varargin)

    [M, N] = size(duomenys);
    Y = duomenys;
    
    [W H] = GetStartingWH(M, N, R);
    
    [W H] = Popt(varargin, 'W0', W, 'H0', H);
    
    ee = 1e-6;
    
    for k=1:niter
        W = W .* (Y*H' )./((W*H)*H' + ee);
        H = H .* (W'*Y )./(W'*W*H + ee);

    end
    
    
end
